# 2024/10/4

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants and parameters (fixed effects for PK and PD models)
Ka = 4.0488   # absorption rate constant
V = 2.561     # volume of distribution
Ke = 0.1439   # elimination rate constant
Tlag = 0.8454 # lag time for absorption

# PD parameters
lambda0 = 0.05   # first-order rate constant of tumor growth
lambda1 = 3.0    # zero-order rate constant of tumor growth
E0 = 340.0       # initial tumor volume
k1 = 0.04        # transit rate constant
k2 = 0.4         # decay rate constant
psi = 20         # switch parameter for inhibition

# Define the PK model
def pk_model(t, y, Ka, Ke):
    Aa, Al = y
    dAa_dt = -Ka * Aa  # Absorption compartment equation
    dAl_dt = Ka * Aa - Ke * Al  # Central compartment equation
    return [dAa_dt, dAl_dt]

# Define the PD model
def pd_model(t, y, C, lambda0, lambda1, k1, k2, psi):
    E1, E2, E3, E4 = y
    E = E1 + E2 + E3 + E4  # Total tumor volume
    Inh = (1 + (lambda0 * E / lambda1)**psi)**(1 / psi)  # Inhibition function
    dE1_dt = (lambda0 * E1 / Inh) - k2 * C * E1  # Proliferating cells
    dE2_dt = k2 * C * E1 - k1 * E2  # Damaged cells (first compartment)
    dE3_dt = k1 * (E2 - E3)  # Transit compartment 2
    dE4_dt = k1 * (E3 - E4)  # Transit compartment 3
    return [dE1_dt, dE2_dt, dE3_dt, dE4_dt]

# Initial conditions  肿瘤初始状态
AaDose = 1.0  # Initial dose
Aa0 = AaDose  # Initial amount of drug in absorption compartment
Al0 = 0.1     # Initial amount of drug in central compartment
E1_0 = 0.7 * E0     # Initial tumor volume in E1
E2_0 = 0.2 * E0  # Initial tumor volume in E2
E3_0 = 0.1 * E0  # Initial tumor volume in E3
E4_0 = 0.0       # Initial tumor volume in E4

# Time points for simulation
t_span = (0, 100)  # Simulation from time 0 to 100 (arbitrary units)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the PK model
pk_sol = solve_ivp(pk_model, t_span, [Aa0, Al0], args=(Ka, Ke), t_eval=t_eval)

# Drug concentration over time (C = Al / V)
C_t = pk_sol.y[1] / V

# Solve the PD model using the concentration from PK model
def combined_model(t, y):
    # Interpolate C from the PK model solution at time t
    C = np.interp(t, pk_sol.t, C_t)
    return pd_model(t, y, C, lambda0, lambda1, k1, k2, psi)

# Initial conditions for PD model
pd_sol = solve_ivp(combined_model, t_span, [E1_0, E2_0, E3_0, E4_0], t_eval=t_eval)

# Extract results
E1_t, E2_t, E3_t, E4_t = pd_sol.y
total_tumor_volume = E1_t + E2_t + E3_t + E4_t

# Plot results
plt.figure(figsize=(12, 6))

# Plot PK results
plt.subplot(1, 2, 1)
plt.plot(pk_sol.t, C_t, label="Concentration (C)")
plt.xlabel("Time")
plt.ylabel("Concentration (C)")
plt.title("PK Model")
plt.legend()

# Plot PD results
plt.subplot(1, 2, 2)
plt.plot(pd_sol.t, total_tumor_volume, label="Total Tumor Volume (E)")
plt.xlabel("Time")
plt.ylabel("Tumor Volume")
plt.title("PD Model")
plt.legend()

plt.tight_layout()
plt.show()